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Abstract 


Atmospheric radiation in the infrared (IR) 8-13 pm spectral region contains a 
wealth of information that is very useful for the retrieval of ice cloud properties from 
aircraft or space-borne measurements. To provide the scattering and absorption properties 
of nonspherical ice crystals that are fundamental to the IR retrieval implementation, we 
use the finite-difference time domain (FDTD) method to solve for the extinction 
efficiency, single-scattering albedo, and the asymmetry parameter of the phase function 
for ice crystals smaller than 40 pm. For particles larger than this size, the improved 
geometric optics method (IGOM) can be employed to calculate the asymmetry parameter 
with an acceptable accuracy, provided that we properly account for the inhomogeneity of 
the refracted wave due to strong absorption inside the ice particle. A combination of the 
results computed from the two methods provides the asymmetry parameter for the entire 
practical range of particle sizes between 1 pm and 10000 pm over wavelengths ranging 
from 8 pm to 13 pm. For the extinction and absorption efficiency calculations, several 
methods including the IGOM, Mie solution for equivalent spheres (MSFES), and the 
anomalous diffraction theory (ADT) can lead to a substantial discontinuity in comparison 
with the FDTD solutions for particle sizes on the order of 40 pm. To overcome this 
difficulty, we have developed a novel approach called the stretched scattering potential 
method (SSPM). For the IR 8-13 pm spectral region, we show that SSPM is a more 
accurate approximation than ADT, MSFES, and IGOM. The SSPM solution can be 
further refined numerically. Through a combination of the FDTD and SSPM, we have 
computed the extinction and absorption efficiency for hexagonal ice crystals with sizes 
ranging from 1 to 10000 pm at 12 wavelengths between 8 and 13 pm. 
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Calculations of the cirrus bulk scattering and absorption properties are performed 
for 30 size distributions obtained from various field campaigns for midlatitude and 
tropical cirrus cloud systems. Parameterization of these bulk scattering properties is 
carried out by using second-order polynomial functions for the extinction efficiency and 
the single-scattering albedo and the power law expression for the asymmetry parameter. 
We note that the volume-normalized extinction coefficient can be separated into two 
parts: one is inversely proportional to effective size and is independent of wavelength, 
and the other is the wavelength-dependent effective extinction efficiency. Unlike 
conventional parameterization efforts, the present parameterization scheme is more 
accurate because only the latter part of the volume-normalized extinction coefficient is 
approximated in terms of an analytical expression. After averaging over size distribution, 
the single-scattering albedo is shown to decrease with an increase in effective size for 
wavelengths shorter than 10.0 pm whereas the opposite behavior is observed for longer 
wavelengths. The variation of the asymmetry parameter as a function of effective size is 
substantial when the effective size is smaller than 50 pm. For effective sizes larger than 
100 pm, the asymmetry parameter approaches its asymptotic value. The results derived in 
this study can be useful to remote sensing applications involving IR window bands under 
cirrus cloud conditions. 
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1. Introduction 


Cirrus clouds form a unique component of the atmosphere and significantly regulate 
the Earth’s radiation energy budget [1-3] because of their large spatial coverage and 
temporal persistence. And yet, the radiative forcing of these clouds is far from being well 
understood because of the wide range of observed cirrus cloud properties such as height, 
optical thickness, particle size, and habit. With the recognition of the importance of cirrus 
clouds, a number of programs have been established to better investigate global cirrus 
properties. For instance, intensive field observations regarding cirrus clouds have been 
conducted as a major component of the First ISCCP Regional Experiment (FIRE) (Phase- 
I in October 1986 and in Phase-II in November-December 1991) [4], the European 
experiments on cirrus (ICE/EUCREX) carried out in 1989 [5]. Additional cirrus field 
campaigns are being planned, such as the Cirrus Regional Study of Tropical Anvils and 
Cirrus Layers (CRYSTAL). Besides these field campaigns, new instrumentation has been 
developed for use in space- or aircraft-borne measurement programs. Such programs 
include the Lidar In-space Technology Experiment (LITE) [6], the MODerate resolution 
Imaging Spectrometer (MODIS) [7], and the future Pathfinder Instruments for Cloud and 
Aerosol Spaceborne Observations-Climatologie Etendue des Nuages et des Aerosols 
(PICASSO-CENA) [8]. The data obtained from these instruments will enhance our 
capability of the detection and retrieval of cirrus clouds on a global scale with a 
substantial temporal span. 

A number of studies have proposed methods to infer global cirrus properties from 
infrared (IR) radiance measurements. One benefit to using IR methods to infer cloud 
properties rather than those that incorporate visible or near-infrared bands is that cloud 
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properties are more consistent between daytime and nighttime conditions. The IR 
techniques are also independent of the sun glint associated with the specular reflection 
over water which is prevalent in daytime data. The IR C0 2 slicing technique is effective 
in inferring cirrus height and effective cloud amount (emittance multiplied by cloud 
fraction) [9,10], Various IR algorithms to infer cirrus particle size have been suggested 
(see, e.g., Ref.[l 1] and [12]), whereas the methods to infer cloud thermodynamic phase 
using data from the 8.52, 11, and 12 fim bands have been demonstrated by Strabala et al. 
[13] and Baum et al. [14]. There have been advances in recent years in the development 
and use of well-calibrated interferometers, such as the High Spectral Resolution Infrared 
Spectrometer (HIS) [15-17]. In coming years, the spectral resolution of IR data measured 
from space will increase greatly when interferometers are launched, such as Infrared 
Atmospheric Sounder Interferometer (IASI) on EUMETSAT and the Geostationary 
Imaging Fourier Transform Spectrometer (GIFTS) and the NPOESS (the National Polar 
orbiting Operational Environmental Satellite System) Aircraft Sounding Testbed- 
Interferometer (NAST-I) [18]. With an interferometer, cloud measurements may be 
recorded at thousands of wavenumbers simultaneously. To date, however, few studies 
have explored the use of infrared interferometer data for cirrus cloud property retrieval. 
One difficulty in the retrieval of cirrus cloud properties from passive airborne- or 
satellite-based radiometric data arises because of the difficulty in determining accurately 
the fundamental scattering and absorption cross sections for ice crystals over a realistic 
range of crystal sizes and shapes at IR bands where the applicability of the ray-tracing 
method breaks down [19]. 
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Cirrus clouds are composed of almost exclusively nonspherical ice crystals, as is 
evident from the observations based on aircraft-borne two dimensional optical cloud 
probes (2D-C) and balloon-borne replicator images (see, e.g., Ref.[20] and [21]). It has 
been shown that the spherical approximation for nonspherical ice crystals in terms of 
equivalent volume or projected-area is inadequate and often misleading, as is illustrated 
by Liou et al. [22], In practice, an incorrect specification of the ice crystal model in 
retrieving the optical thickness of cirrus clouds from satellite-borne measurements can 
lead to an underestimation or overestimation of the actual optical thickness by a factor 
that can exceed 3 [23]. Thus, nonsphericity of ice crystals must be accounted for in the 
development of a reliable retrieval algorithm. On this specific issue, the significance of 
using reliable single-scattering properties of ice crystals to generate look-up tables for 
retrieval implementation and the parameterization of the bulk radiative properties of 
cirrus clouds has been demonstrated and articulated in a number of recent publications 
(e.g., Ref.[14, 24-31] ). Although substantial advancements have been made in the 
fundamental study regarding scattering and absorption by ice crystals, as recently 
reviewed by Liou and Takano [24], Mishchenko et al. [34], and references cited therein, 
there is no a single method that can cover the entire size parameter spectrum for light 
scattering computations. To derive the scattering and absorption properties of 
nonspherical ice crystals, Liou et al. [22] developed the concept of a unified theory. This 
unified theory is based on a combination of a numerically accurate finite-difference time 
domain technique (FDTD) for small particles [35-40] and an improved geometric optics 
method (IGOM) [41-42] for large particles at visible and near-infrared wavelengths. 
However, for IR wavelengths where strong absorption is involved, there is a discontinuity 
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between the FDTD and IGOM results at size parameters on the order of 20, which is in 
practice the computational limit for the FDTD method given current computer resources. 
One cause of this discontinuity is the tunneling or the above -edge effect [43]. Several 
approximate methods [29,30, 44, 45] suggested to account for the tunneling effect are 
based on parameterizations involving Mie theory, a combination of Mie theory and the 
geometric optics method, or the complex angular momentum theory developed by 
Nussenzveig and Wiscombe [46,47]. 

The intent of this study is to develop methods capable of deriving the fundamental 
scattering and absorption cross sections for nonspherical ice crystals spanning a size 
range of 1 pm to 10000 pm at IR wavelengths ranging from 8 pm to 13pm, where 
absorption effects can be appreciable. In section 2, we present methodology for 
calculating the scattering properties at infrared wavelengths. In particular, we present the 
stretched scattering potential method (SSPM) to calculate the extinction and absorption 
cross sections. In Section 3, we present the bulk radiative properties and parameterization 
for cirrus clouds. Finally, the conclusions of this study are given in Section 4. 

2. Computation of Optical Properties for Ice Crystals at Infrared Wavelengths 

2. a. Aspect Ratio for Ice Crystals and Selection of Wavelengths 

A variety of nonspherical ice crystal habits, or shapes, have been observed in cirrus 
clouds, including hexagonal columns and plates, bullet rosettes, and complex aggregates. 
To demonstrate our methodology, we simplify the complexity of determining what best 
represents the habit distribution in cirrus by assuming cirrus are composed solely of 
hexagonal shapes with a random orientation in space. Studies of in-situ cirrus 
measurements have demonstrated that the upper layers of midlatitude cirrus cloud 
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systems are often comprised of pristine hexagonal crystals [48]. Hexagonal ice crystals 
are capable of producing 22° and 46° halos and other optical phenomena by scattering 
incident solar radiation at short wavelengths. The hexagonal shape model has been often 
assumed in previous studies concerning cirrus clouds [24-30, 49-51]. 

To carry out scattering calculations, we define the aspect ratio for ice crystals at 
various sizes as follows: 


' 1, L < 40 \im 

2 a!L = < exp[-0.017835(L— 40)], 40 < L < 50 pm 
5.916/L 1/2 ,L>50 pm, 


( 1 ) 


where a is the semi-width of cross section and L is the length of an ice crystal. The aspect 
ratio defined by Eq.(l) roughly corresponds to the observations reported by Ono [52] and 
Auer and Veal [53]. Figure 1 shows a comparison of the aspect ratio defined in this study 
with that of Takano and Liou [49], who used only five size-bins in discretizing the size 
distribution for ice crystals. Excellent agreement is evident for the comparison. In the D- 
L plot of Fig. 1 in which D=2a is the width of cross section, the continuous aspect ratio is 
not smooth at a size of L = 40 pm because of an abrupt variation in the derivative of 
aspect ratio. However, this effect is very small. 

To economize the computational effort, the scattering calculations are performed 
for 12 wavelengths between 8 pm and 13 pm, which are selected at 8.0, 8.5, 9.0, 9.5, 
10.0, 10.5, 10.8, 11.0, 11.5, 12.0, 12.5, and 13.0 pm. The choice of these wavelengths is 
based on the characteristics of ice refractive index as shown in Figure 2 [54], The dotted 
vertical lines indicate the locations of the selected wavelengths. Note that the variation of 
the refractive index is essentially linear in the intervals of the selected wavelengths. 
Because the refractive index is essentially linear in a given interval, the scattering and 
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absorption properties at any arbitrary wavelength in the 8-13 pm region can be 
approximated by an interpolation of those calculated at the selected wavelengths. We also 
note that near 10.8 pm there exists a region of extremely strong absorption, often called 
the Christiansen band [55,56]. Near 10.8 pm, the real part of the refractive index has a 
value close to 1 and the extinction efficiency reaches a minimum. As is evident from 
Fig.2, the imaginary part of the refractive index is significant at the Christiansen band, 
leading to the dominance of the absorption effect in the extinction process. The behavior 
of ice absorption in the Christiansen band has been shown in laboratory measurements 
[56,57] and a theoretical explanation was reported by Yang et al. [58]. 

2.b. Asymmetry parameter for angular distribution of scattered energy 

The angular distribution of the energy scattered by a particle is defined by its 
scattering phase function. At IR wavelengths, the phase function of an ice crystal is 
essentially featureless when compared to that at a visible or near-infrared wavelength. 
The primary scattering features of a hexagonal crystal at a visible wavelength, such as the 
halo peaks and backscattering enhancement, are smoothed out at IR wavelengths because 
of absorption. 

For radiative transfer calculations at IR wavelengths, the detailed information of the 
phase function is unnecessary because the multiple scattering effect is only on the order 
of a few percent. Thus, the asymmetry parameter that describes the magnitude of the 
deviation of particle phase function from isotopic scattering is very useful at IR 
wavelengths, which is defined as follows: 

g = — f P(0)cos0sin0</0, (2) 

2 J 0 
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where P(0) is the normalized phase function. Given the asymmetry factor associated 
with an IR phase function of nonspherical ice crystals, an accurate phase function for a 
given size distribution of ice crystals can be approximated by the Henyey-Greenstein 
analytical function: 

P " G(0) = a + * 2 -2*cos0) 3 ' 2 

N 

»^(2/ + l)g / /}(cose), (3) 

1=0 

where/} (cos 0) is the set of Legendre polynomials derived from decomposition of the 
phase function. The advantage of using the Henyey-Greenstein phase function in 
radiative transfer calculations lies in its simplicity and efficiency in the expansion of a 
phase function in terms of Legendre function. There is a physical justification for such a 
simplification. For an IR wavelength, the phase function calculated for a given size 
distribution of ice crystals tends to be smooth at sidescattering and backscattering angles 
whereas a strong diffraction peak is typically noted in the forward scattering direction. 
The forward peak may need to truncated; for this the delta-M method developed by 
Wiscombe [59] may be employed. The truncated phase function can be well 
approximated by the analytical Henyey-Greenstein function at the IR wavelengths 
relevant to our study. 

The FDTD technique is employed to solve for the scattering and absorption 
properties for ice crystals whose maximum dimensions are smaller than 40 pm. The 
FDTD method solves the interaction of electromagnetic waves with a particle of any 
given shape and is based on the difference analog of time-dependent Maxwell equations. 
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The strength of this method is its simplicity in concept, computational flexibility, and 
robust nature in dealing with light scattering problems involving nonspherical and/or 
inhomogeneous particles. The disadvantage of this method is its tremendous demand on 
computational resources when the size parameter is larger than 20. Detailed descriptions 
of the FDTD method may be found in Yang and Liou [38]. 

To calculate the asymmetry parameter for ice crystals larger than 40 pm, we use a 
hybrid method based on the improved geometric optics principle (IGOM) and the 
electromagnetic integral equations [41, 42]. At IR wavelengths, the refractive wavelength 
inside an ice crystal is inhomogeneous because the planes of constant phase are not 
parallel to the planes of constant amplitude [60]. This inhomogeneity will affect both the 
ray direction and the reflection/refraction components at the air-particle interface [35]. In 
particular, the Fresnel coefficients are not unique when absorption is involved, as was 
shown recently by Yang et al. [61], who further determined the proper form of the 
refraction/refraction coefficients that should be used in ray-tracing calculations. We have 
accounted for this improvement in the present set of scattering calculations. 

Figure 3 shows the asymmetry parameter calculated by the FDTD and the IGOM 
over a range of ice crystal sizes at wavelengths 8.5, 11, and 12 pm. The radiative 
information at these three wavelengths are being used to retrieve cloud thermodynamic 
phase in MODIS data. For comparison of results between spherical and nonshperical 
particles, results computed from Mie theory for equivalent ice spheres are provided in 
Fig. 3. Following work by Mitchell and Arnott [62], Fu et al. [29], Yang et al. [63], and 
Grenfell and Warren [64], we define the radius of the equivalent sphere for an individual 
hexagonal ice crystal as follows: 
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( 4 ) 


_ 3 volume _ V 3a 2 L 

e 4 projected- area V3a 2 / 2 + aL ’ 

where a and L the semi-width and the length of ice crystal, respectively. The advantage 

of the preceeding definition over those for equivalent volume (or projected-area) is that 

both volume and area are conserved in approximating the nonspherical particle in terms 

of a sphere. The inserted sub-diagrams in Fig.3 provide enlargements for ice crystal sizes 

of 1-40 pm. For small ice crystals with sizes of 3-9 pm, the asymmetry parameters 

computed from Mie theory for the equivalent spheres and the FDTD technique for the 

nonspherical particles are quite different. For sizes larger than 100 pm, the IGOM results 

essentially converge to the Mie results. This convergence occurs because the rays 

refracted into the particles are almost entirely absorbed and the diffraction and external 

reflection contributions dominate the scattered energy. Under random orientation 

conditions, diffraction and external reflection are not sensitive to the details of particle 

shape. 

From Fig.3, we note that the FDTD results and IGOM solutions converge at a particle 
size of 40 pm. Although the absolute amount of scatted and absorbed energy calculated 
by IGOM can have a substantial error for moderate particle sizes on the order of 40 pm 
or less, the IGOM can be used to predict the relative angular distribution of scattered 
energy. Thus, the normalized IGOM phase function can be a good approximation. From 
this physical rationale and the numerical results shown in Fig.3, the concept of the unified 
theory developed by Liou et al. [22] can also be applied to the computation of the 
asymmetry parameter. Therefore, we combine the FDTD solution (for sizes smaller than 
40 pm) and IGOM results (for sizes larger than 40 pm) for the asymmetry factor so that 
the entire size spectrum can be covered. 
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Figure 4 shows merged FDTD/IGOM results for the asymmetry factor for the 12 
wavelengths for particle sizes of 1-10000 pm. The variation of the asymmetry factor 
with the particle size displays a minimum at 30 pm for wavelengths 8.0, 8.5 and 9.0 pm 
where the absorption of ice is moderate. This minimum vanishes for longer wavelengths. 
The overall feature of the asymmetry parameter is that the g values are close to unity for 
large sizes because diffraction dominates the scattered energy which is concentrated in 
forward direction. For all 12 wavelengths, there is a substantial increase in the value of 
the asymmetry parameter as the size of the particle increases from 1 pm to 40 pm. For 
sizes larger than 40 pm, the asymmetry parameter essentially reaches an asymptotic 
value. 

2.c. Stretched Scattering Potential Method for Computing Extinction and Absorption 
Cross Sections for Nonspherical Particles 

Several methods exist to calculate extinction and absorption cross sections for 
nonspherical particles. To date, the upper limit of the size parameter region for which one 
can obtain the exact solution for the scattering properties of certain nonspherical particles 
is on the order of 200. This limitation is for axisymmetrical particles such as spheroids, 
finite circular cylinders, and so-called Chebyshev particles [65], as solved by T-matrix 
method [66]. The applicable size parameter regimes for other exact methods to solve the 
scattering properties of nonspherical particles are normally smaller than that associated 
with the T-matrix method. For a non-axisymmetrical particle, the size parameter region 
for which an exact solution may be obtained is substantially reduced. For instance, the 
limitation on size parameter for the T-matrix method is on the order of 40 when it is 
applied to hexagonal particles [A. Baran, personal communication]. The FDTD and its 
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counterpart, the discrete dipole approximation (DDA) [67,68], encounter computational 
difficulties for size parameters larger than 20. For these reasons, use of an approximate 
method is necessary in practice. The anomalous diffraction theory (ADT) developed by 
van de Hulst [69] has been widely used to calculate the extinction and absorption cross 
sections for nonspherical particles [70-72]. Unfortunately, ADT leads to substantial errors 
if the refractive index of the particles is not close to 1. In particular, ADT fails to account 
for the tunneling effect; the absorption efficiency calculated from ADT cannot be larger 
than 1. The intent of this section is to derive an improved approximate method to 
calculate the extinction and absorption cross sections. 

The wave equation derived from the Maxwell equations can be written as 

(y 2 + k 2 )E(r) = -U(r)E(r), (5) 

where U(r) = [m 2 (r)-l]k 2 in which m is the complex refractive index of the medium 
and k is the wavenumber in a vacuum. Through comparison of Eq.(5) to the standard 
Schrodinger equation in quantum mechanics, we can regard U(r ) as a scattering 
potential. Thus, the light scattering problem can be treated as one concerning the 
interaction of photons and a potential. A separation of variables method can be used to 
solve Eq.(5) exactly but is limited to only a few particle shapes such as spheres and 
spheroids. In practice, a high energy approximation (HEA) [73-75] or ADT can be used 
to solve Eq.(5) approximately. These methods are the eikonal type [76] in which 
rectilinear projectiles are assumed for the propagation of photons, with (or without) 
deviation only at the places where the scattering potential vanishes, i.e., at the particle 
surface. Thus, in these conventional approximations, the nonzero-interaction region 
where wave function undergoes phase delay and absorption is limited inside the particles. 
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Recently, Fu et al. [30] have computed the Poynting vector for the near-field associated 
with the scattering of IR radiation by a hexagonal ice crystal to illustrate the flow of 
electromagnetic energy around the particle. It is shown that the incident Poynting stream 
outside the geometric projected area of the particle can deviate and be traced through the 
particles, a phenomenon otherwise known as the tunneling or the above-edge effect. 
Although HEA and ADT provide a good physical insight in concept and a simple 
mathematical formulation in practice, they may produce significant errors in 
computation, e.g., a substantial underestimation of absorption efficiency in the resonance 
region because of the failure to account for the tunneling effect. 

To avoid this shortcoming of the conventional eikonal-type approximation that is 
applied to solve Eq.(5), it is necessary to stretch the scattering potential to account for 
the tunneling effect. Figure 5(a) illustrates the region of non-null scattering potential, 
where the interaction of photons and medium occurs for the conventional eikonal 
methods. The scattering potential U(r ) is nonzero only inside the particles for the 
conventional methods. Figure 5(b) is the conceptual diagram for the present stretched 
scattering potential method (SSPM), where a is the physical size of the particle and Aa is 
extension for the stretched potential. In this method, the scattering potential U{r) is 
stretched so that the non-zero effect region for photon propagation extends outside the 
physical volume of the scattering particle. We assume the scattering potential has a 
quadratic distribution, for example, in the spherical particle case, as follows: 

(m 2 -l)k 2 , r<a 

U(r) = j (m 2 — l)fc 2 [l - (r-a)/ Aa] 2 , a<r<a + Aa. (6) 

0, r < a + Aa 
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For stretched scattering potential in the region between a and a+Aa, we can define an 


effective refractive index m = ( m r + im [ ) given by 

m 2 -l = (m 2 -l)k 2 [l-(r-a)/Aa] 2 . (7) 

The solution of Eq.(7) is 

a = 1 + (m 2 -mf - 1)[1 -(r-a)/ Aa] 2 , (8a) 

P = m r m l [1 -(r-a)/ Aa] 2 . (8b) 

m r (r) = [a + (a 2 +4p 2 ) 1/2 ]/2, and (8c) 

m i {r) = $/m r {r), (8d) 


The external reflection and multiple internal reflections and subsequent transmission of 
the wave is completely neglected in the conventional eikonal-type approximate methods 
to solve Eq.(5). The result of ignoring these effects is equivalent to an overestimation of 
scattering potential. To include the external and internal reflection/transmission effect 
equivalently, the scattering potential should be compressed instead of being stretched. 
The compressed scattering potential is given in the form of 

(m 2 -l)k 2 , r<a-Aa 

U(r) = < (m 2 -l)k 2 [l-(r- a + Aa)/ Aa] 2 , a- Aa<r<a. (9) 

0, r < a 

Figure 6 shows the concept of using an eikonal approximation in the ADT framework 
to calculate the extinction and absorption cross sections for a spherical particle. In 
practice, continuous variation of the effective refractive index is discretized as various 
layers. That is, the effect of stretching the scattering potential is equivalent to adding 
some dielectric medium layers outside the scattering particle. For a given layer, we 
calculate the mean effective refractive index as follows: 
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( 10 ) 


1 /•r > +Ar/2 

<m r + im i > j— — [m r (r) + im^r^dr, 

Jr r &rl2 

where r y is the radius of the ;'th layer and A r is the thickness of the layer. A uniform 

thickness is assumed for all the layers in the present study. According to the geometry 
shown in Fig. 6 and the physical assumption in conjunction with Huygen’s principle for 
the ADT approximation [69], we obtain expressions for the extinction cross section 
( <5 ext ) and the absorption cross section (d abs ) for a layered sphere, respectively, as 

follows: 

o ext = 2jj(l-e~ y cosp )d 2 P, (11a) 

p 

Cabs=\\d-e~ 2y )d 2 P, (lib) 

p 

where p is the phase delay that a photon undergoes in conjunction with its penetration 
of the multilayered sphere. The penetration parameter is given by 

N 

p = kd Q (m r -\) + k(dj - dj _\ )(< m r >y - 1) , (11c) 

i=l 

in which N is the number of total layers outside the particle. The attenuation factor, y , for 

the damping of the incident wave is given by 

N 

y = kdQm i + k(dj - dj_ x ) < in t >j. (lid) 

M 

Eqs.(l lc) and (1 Id) are for the case where the scattering potential is stretched outside the 
particle. A similar mathematical formulation may be derived for the case when the 
scattering potential is compressed or reduced. 
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From the preceding discussion of the SSPM approximation, the only tuning 
parameter is the distance Aa for stretching the scattering potential. The physical 
processes associated with the extinction and absorption of the incident wave are quite 
different. For example, the surface wave (a term in the rigorous physical picture) [43, 69] 
that creeps along the particle surface can contribute to the extinction but not to the 
absorption. Therefore, the magnitude of the extension of the scattering potential may be 
different for calculating the extinction cross section than for the absorption cross section. 
For this reason, we denote Aa as A a e and A a a for extinction and absorption calculations, 

respectively. 

Figure 7 shows the values of A a e that are derived from the best fit of the exact 

Mie solution using the SSPM results in the spherical case. In this study, the 
computational Mie code developed by Wiscombe [77] is used. To determine A a e , a Mie 

calculation is first carried out; subsequently a Monte Carlo method is employed to 
determine the proper value for A a e for the SSPM method. Specifically, for each SSPM 

calculation in the case where the potential is stretched rather than compressed, we chose 

Aa e =\a, (12) 

where I; is a random number that is uniformly distributed in [0,1]. From Fig.7, it is shown 
that the magnitude of A a e increases with particle dimension. However, the ratio of A a e 
to the particle dimension reaches a maximum at a particle size of approximately 10 pm 
where the resonance effect is largest. As is evident in Fig.7, the SSPM extinction 
efficiency actually overlaps with the Mie solution. If the procedure for calculating A a e 

and the extinction efficiency were reversible, the SSPM could reproduce the exact theory, 
providing that adequate extension of the scattering potential was known a priori. It is 
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straightforward to obtain the correct stretching scattering potential for the SSPM 
calculation for a spherical particle case because of the availability of the exact Mie 
theory, as is evident from the results shown in Fig.7. For scattering by a nonspherical 
particle, it is necessary to define an approximation of the correct expansion of the 
scattering potential. 

Figure 8 is similar to Fig. 7 except it relates to the absorption calculation. Again, 
if A a a is properly selected, the SSPM result can match the analytical Mie solution. 
Another feature of the behavior of A a a as a function of the particle size is that negative 
values are noted for very large particles. This means that the scattering potential should 
be compressed in the calculation of absorption cross section because the effect of external 
reflection and refraction is not accounted for in the eikonal type approximation given by 
Eqs.(lla)-(lld). For the absorption efficiency at 12.0 pm wavelength, a pronounced 
tunneling effect can be observed because the extinction efficiency is substantially larger 
than unity for particle sizes near 10 pm. The SSPM accounts for the tunneling effect, 
providing a proper expansion of the scattering potential is used. 

It is problematic to accurately predescribe the expansion of the scattering 
potential for SSPM calculation for a nonspherical particle. In this case, we approximate 
A a e and A a a by the values obtained from the equivalent spheres with radii defined by 

Eq.(4) that conserves both the volume and projected-area in the equivalence process. In 
the SSPM computation for hexagonal ice crystals, we specifically solve for the phase 
delay and wave attenuation for a number of layered hexagons. The total number of the 
hexagonal layers and their thicknesses are approximated by their counterparts in the 
spherical case. Figure 9 shows the absorption efficiency of hexagonal ice crystals 
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calculated from various methods. The errors of the FDTD results are less than 1% for the 


typical grid resolution [40]. Thus, the FDTD results here can be used as a reference for 
checking the accuracy of the approximate methods. As shown in Fig. 9, the SSPM 
results essentially overlap with the FDTD solution for wavelengths of 8.0, 8.5, 9.0, 9.5, 
10.0, and 10.5 pm. For longer wavelengths, deviation of the SSPM results from FDTD is 
noted for particle sizes larger than 20 pm. It also becomes evident that the GOM, ADT, 
and MTFES may lead to substantial errors. In particular, ADT will significantly 
overestimate the absorption efficiency for very large particles when the particles are 
strongly absorptive at a given wavelength. The asymptotic value for ADT absorption 
efficiency is unity because the external reflection is unaccounted for in this 
approximation. The Mie results shown in Fig. 9 are generated for spheres whose radii are 
defined by Eq.(4). 

Figure 10 is similar to Fig. 9 except that results are shown for extinction 
efficiency. One obvious result is that Mie theory for the equivalent spheres substantially 
overestimates the extinction at the resonance maximum. The GOM and ADT methods 
converge for large particles. The GOM results shown here are from the improved 
geometric optics method developed by Yang and Liou [41, 42]. If the conventional ray- 
tracing method is used, the extinction efficiency is simply 2 regardless of particle size, as 
is pointed out by Yang and Liou [35]. For all 12 wavelengths the SSPM results seem to 
provide the most consistent comparisons with the FDTD method. 

Some physical processes, such as the external reflection and detailed nonsphericity 
effects, are not fully accounted for in the SSPM because exact expansion of the potential 
is not used. The effects associated with these physical processes are reflected in the 
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equivalent Mie and GOM results to some extent. Thus, we can refine the SSPM results as 
follows: 

Qrefined = 0 “ 5 2 ) QssPM + $1 @Mie + §iQgOM' ( 13 ) 

where Q stands for either extinction efficiency or absorption efficiency Q a . The 
values of the parameters 8] and 5 2 are so determined that the G re f me d best fit the FDTD 

results for crystal sizes between 20 pm and 40 pm. The refined SSPM results are given 
only for sizes larger than 20 pm, because the geometric optics solution involved in 
Eq.(13) is essentially meaningless due to the failure of the localization principle for small 
size parameters. Figure 11 shows the refined SSPM results for the extinction and 
absorption cross section efficiencies. With this refinement, calculations can be provided 
for the entire size spectrum. The refinement procedure employed here is similar to that 
used by Fu et al. [29] who use a weighted summation of Mie and GOM solutions to 
obtain an approximate fit for the single-scattering properties of hexagonal ice crystals. 

3. Application to Cirrus Clouds and Parameterizations 

To derive a set of ice crystal bulk scattering properties for practical applications, it 
is necessary to derive the single-scattering properties using realistic particle size 
distributions. In this study, we select the 28 size distributions used by Fu [28] with an 
additional two datasets from Mitchell et al.[78]. These size distributions were measured 
for a variety of midlatitude and tropical cirrus clouds and were obtained during various 
field campaigns including the Central Equatorial Pacific Experiment (CEPEX) [79]. 
McFarquhar et al. [80] showed that the microphysical properties of ice crystals for 
tropical cirrus clouds are quite different from those for midlatitude cirrus systems. 
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It is common to characterize the bulk properties of a size distribution by two 
parameters - ice water content ( IWC) and effective size. For a given size distribution, 
IWC is defined as 


IWC 


- p ,«J 


V(L)n(L)dL, 


(14) 


where p ice is the mass density of bulk ice, n(L ) is particle number density, V{L) is the 
volume of an ice crystal with maximum dimension of L, and L min and L max are the 
minimum and maximum sizes in the size distributions, respectively. There are many 
definitions for effective size in the literature. Wyser and Yang [31] performed a 
comprehensive comparison of definitions commonly used in different parameterization 
efforts. In this study, we define the effective size following Foot [81], Francis et al. [5], 
Fu [28], Wyser and Yang [31], and Grenfell and Warren [64] as 


D e = 


,C 

i: 


V(L)n{L)dL 


A(L)n(L)dL 


(15) 


where A(L) is the projected-area of the particle with size of L. The preceeding definition 
reduces to that defined by Hansen and Travis [82] for spherical particles. We note that 
Eq.(15) is similar to Eq.(4) except that only an individual particle is considered in Eq.(4) 
whereas Eq.(15) considers a size distribution. In addition to the effective size, we also 
define the mean maximum dimension for a given size distribution as follows: 


<L> 


i 


Ln{L)dL 


^tnin 


J 


(16) 


n{L)dL 
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Figure 12 shows a scatterplot for the effective size D e as a function of the mean 
maximum dimension <L> for the aforementioned thirty size distributions described. It is 
evident that D e is normally larger than <L> for a given size distribution. This occurs 
because the distributions contain a large number of small particles that are significant in 
defining the mean maximum dimension. However, the contribution of small particles to 
the total volume or projected-area is relatively small regardless of their contribution. 
Thus, the largest weight in defining D e is for moderate or large ice crystals. It is also 

evident from Fig. 12 that the spectrum of effective sizes ranges from approximately 10 
pm to 160 pm for the selected size distribution datasets. Given the range of effective 
sizes, the present set of ice cloud size distributions should provide an effective basis for a 
parameterization effort. 

For a given size distribution of ice crystals, the bulk extinction coefficient is 
defined as follows: 


J, 


P e = QtiL) A{L)n{L)dL. 


(17) 


It is convenient in practice to define a volume-normalized extinction coefficient [31], 
given by 


K = 




t 


Q,{L)A(L)n{L)dL 


( IWCIp ice ) 


J 


(18) 


V(L)n(L)dL 


Note that the quantity ( IWC/p ice ) in the preceeding equation provides the total bulk 
volume of ice crystals for a given size distribution. For a cirrus cloud with a given IWC 
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and geometrical thickness A z, the dimensionless optical thickness of the cloud may be 
obtained from 

~ IWC 

T = Az. (19) 

Pice 

This relationship is useful in practice because IWC is now a prognostic parameter in 
many global climate models (GCMs) [83]. In many parameterization efforts, it is 
common to parameterize P e or as a function of the effective size. In particular, it has 
been found that the extinction coefficient decreases with an increase in effective size 
(e.g.,Ref.[28,31] ). We suggest that it may be more useful to parameterize the mean 
extinction efficiency rather than the extinction coefficient. The volume-normalized 

extinction coefficient P g can be rewritten exactly as 

P«=4<a>. (20 

where < Q e > is the mean extinction efficiency that is defined as 

f " Q e (L)A(L)n(L)dL 

<Qe>=- h Jl • (2D 

A(L)n(L)dL 

From inspection of Eqs. (20) and (21), we suggest that the accuracy of a parameterization 
scheme would be improved if the mean extinction efficiency could be parameterized 
accurately. 

For visible wavelengths, the size parameter is very large even for small ice crystals on 
the order of 10 pm. Thus, for visible wavelengths, the extinction efficiency of an ice 
crystal essentially approaches a geometric optics asymptotic value of 2 for a realistic set 
of particle sizes. For the IR region of interest in the present study, the resonance effect is 
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obvious for ice crystal sizes between 10 |im and 20 pm. However, the number 
concentration of particles in this size region may be significant for some ice clouds such 
as contrails and cold cirrus. Thus, <Q e > needs to be investigated carefully in a 

parameterization effort. From Eq.(20), the wavelength-dependent part of the extinction 
coefficient is the mean extinction efficiency. The mean extinction efficiency can be 
parameterized as a function of effective size. 

Figure 13 shows the variation of the volume-normalized extinction coefficient, 
as a function of effective size for 12 wavelengths in the IR window region. The 
extinction coefficient decreases with an increase of effective particle size, as expected 
from Eq.(20). For effective particle sizes larger than 80 pm, the extinction coefficient 
approaches its asymptotic value with a small variation as particle size increases. We note 
that P e can be regarded as the extinction coefficient for a unity IWC because the mass 
density of ice is independent of the size distribution. For a fixed cloud geometrical 
thickness and IWC, according to Fig. 13 and Eq.(19), the cloud will have a larger optical 
thickness if the ice crystals are small. Conversely, the optical thickness of the cloud will 
be reduced if the particle sizes are large. 

The circle symbols in Fig. 14 show the variation of the dimensionless mean 
extinction efficiency as a function of H D e . The variation pattern of <Q e > is shown to 

have a dependence on wavelength. The magnitude of <Q e > ranges between 1.5 and 2.2. 
When 1 1 D e approaches zero (i.e., large effective size), the mean effective extinction 

efficiency is close to 2. This is expected because the particle extinction efficiency 
reaches its geometric optics asymptotic value when the size is very large. For 
wavelengths 8.5, 9.0, 9.5, 10.0, 12.0, 12.5 and 13.0 pm, an extinction maximum is 
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observed for the moderate effective sizes. For wavelengths 10.5, 10.8, 11.0, and 11.5 pm 
near the Christiansen band, the mean extinction efficiency decreases with the increase of 
1 / D e , i.e., a decrease of particle mean size. 

The mean extinction efficiency may be parameterized by a second-order polynomial 
in the form of 

<Qe>= n 2 (l/^) 2 + T li(l / ^) + Tlo> ( 22 ) 

where the coefficients T| 0 , rfo, and r) 2 are determined from the least-squared fitting 
technique. Table 1 lists the three fitting coefficients derived for the 12 wavelengths. The 
parameterization results are also shown in Fig. 14 (solid lines). 

The circle symbols in Fig. 15 show the variation of single-scattering albedo with 
1 !D e . Here the mean value of the single-scattering albedo for a given size distribution is 

defined as follows: 

»L 

I nux 

J [Qe(L)-Q a (L)]A(L)n(L)dL 

<(0>= — ^ — i , (23) 

[ Qe(L)A(L)n(L)dL 

The overall variation trend of < cb > in the IR window region can be grouped into two 
catagories: for wavelengths smaller than 10.0 pm, the single-scattering albedo increases 
with an increase of 1 / D e , i.e., single-scattering albedo increases with decreasing mean 

particle size; for wavelengths longer than 10.5, the opposite behavior is observed. A 
second-order polynomial function may be employed to fit the single-scattering albedo in 
the form of 

<©>=C 2 (1/^) 2 +Ci(1/^) + ^o- (24) 
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Table 2 lists the fitting coefficients in Eq.(24). The solid lines in Fig. 15 show the 
corresponding parameterization results. 

The circle symbols in Fig. 16 show the asymmetry factor calculated for the 30 
size distributions at 12 wavelengths. In this figure, the asymmetry factor is plotted against 
the effective size rather than against 1 1 D e as in Figs. 14 and 15. The asymmetry factor 
increases with an increase in the effective size. One explanation for this is that particle 
absorption reduces the amount of scattered energy in the side scattering and back 
scattering directions, and affects the transmission of the incident wave. In addition, the 
diffracted energy is concentrated in a narrower region around the forward scattering 
direction when the particle sizes are increased. 

While polynomial functions are employed often to parameterize the asymmetry 
factor (e.g.,Ref.[28,31]), we suggest instead the use of a power law form, which may be 
mathematically expressed as follows: 

< g>= £(< D e >) K * < 25 ) 

The fitting coefficients <; and K for Eq.(25) are listed in Table 3. Our parameterization 
results as given by Eq.(22), (24) and (25) should be limited to the domain of effective 
sizes (approximately 10-160 pm) used in our analyses. For an effective size outside of 
this range, the applicability of the parameterization should be checked discreetly to avoid 
any unpleasant artifacts. 

Conclusions 

We present fundamental scattering and absorption properties for hexagonal ice 
crystals with sizes ranging from 1 pm to 10000 pm in the infrared (IR) 8-13 pm window 
region. The 8-13 pm region contains a wealth of spectral information that may be 
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exploited by satellite-borne IR retrieval of ice cloud properties. To obtain scattering and 
absorption properties over this range of particle sizes, two methods need to be employed 
because there is currently no single scattering computational method that can cover the 
entire size parameter spectrum. In this study we use the finite-difference time domain 
method to solve for the extinction efficiency, single-scattering albedo, and the asymmetry 
parameter of the phase function for ice crystals smaller than 40 pm. We find that the 
improved geometric optics method (IGOM) can be employed to calculate the asymmetry 
parameter for ice crystals larger than 40 pm if one accounts for the inhomogeneity effect 
of the refracted wave inside the nonspherical ice particle. The combination of the two 
methods provide the results for the entire range of particle sizes over wavelengths 
ranging from 8 pm to 13 pm. 

Other methods are compared but shortcomings are noted. In particular, the 
geometric optics and the anomalous diffraction theory fail to account for the tunneling 
effect (a phenomenon that the incident energy stream outside the particle projected-area 
can be trapped and scattered/absorbed) in calculating the absorption and extinction 
efficiency. The analytical Mie theory does account for the tunneling effect. When the 
equivalent spherical approximation is applied to scattering by a nonspherical ice crystal, a 
significant overestimation of the effect can result. This is a particular problem at the size 
parameter where the resonance maximum is produced. For the extinction and absorption 
efficiency calculations, several methods including the IGOM, the Mie solution for 
equivalent spheres, and the anomalous diffraction theory (ADT) can lead to a substantial 
discontinuity at a particle size of 40 pm. 
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We have developed a novel new approach called the stretched scattering potential 
method (SSPM) to address the aforementioned difficulties. At 12 wavelengths in the 
spectral region 8-13 pm, we show that the SSPM is a more accurate approximation than 
ADT, Mie theory, and IGOM. Additionally, we suggest further numerical refinements to 
the SSPM solution. Through a combination of the FDTD and SSPM, we have computed 
the extinction and absorption efficiency for hexagonal ice crystals with sizes ranging 
from 1 pm to 10000 pm at 12 wavelengths between 8-13 pm. 

Based on the single-scattering properties obtained for individual ice crystals, 30 
size distributions obtained from various field campaigns for midlatitude and tropical 
cirrus cloud systems have been selected to calculate the bilk scattering/absorption 
properties for the clouds. A further parameterization effort is carried out to analytically fit 
these bulk scattering properties by using second-order polynomial functions for the 
extinction efficiency and the single-scattering albedo and a power law form for the 
asymmetry parameter. We note that the volume-normalized extinction coefficient can be 
separated into two parts: one is inversely proportional to effective size and is independent 
of wavelength, and the other is mean effective extinction efficiency. Unlike conventional 
parameterization efforts, in the present parameterization scheme the latter part of the 
volume-normalized extinction coefficient is numerically fitted. Finally, it should be 
pointed out that single-scattering properties of ice crystals in the IR window region show 
a very strong wavelength-dependence, in particular, in the case of single-scattering 
albedo. Numerical results show that the single-scattering albedo increase with the 
decrease of particle sizes for wavelengths shorter than 10.0 pm and opposite pattern is 
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observed for longer wavelengths. The results obtained in this study can be useful for 
spaceborne IR retrieval of cirrus clouds. 
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Figure Captions 

Fig. 1 The aspect ratio used for defining the geometry of the hexagonal ice crystals. D 
and L are the width of cross section and length of the particle, respectively. 

Fig.2 The variation of refractive index in IR (7.5-13.5 pm) region. The refractive 
index data are those compiled by Warren [1984]. The dotted lines in the 
diagrams indicate the locations of the wavelengths selected for the presented 
scattering computation. 

Fig.3 The asymmetry parameter computed from the FDTD technique and the 
improved geometric optics method for hexagonal ice crystals at three 
wavelengths. The subdiagrams are the enlargements for particle sizes between 1 
to 40 pm. The Mie solutions are for the equivalent spheres whose radii are 
defined by Eq. (4). 

Fig.4 The combined FDTD and geometric optics solutions for the asymmetry factor at 
12 wavelngths, which cover the particle sizes from 1 to 10000 pm. 

Fig.5 The conceptual diagrams in spherical case to illustrate the principle of the 
conventional eikonal-type approximation and the present stretching scattering 
potential method (SSPM). (a): the non-zero effect region for phase delay and 
wave attenuation in the conventional method. The radius of the sphere is 
denoted by a. (b): the region of nonzero scattering potential in the SSPM. The 
scattering potential is expanded outside the particle in the region between a and 
a + Aa to account for the tunneling effect. 

Fig.6 The conceptual diagram to illustrate the computational scheme for multi-layered 
particle in the frame of an eikonal-type approximate method. Rectilinear 
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projectiles are assumed for the of photon paths in conjunction with the 
computation of phase delay and wave attenuation. 

Fig.7 The expansion parameter A a e derived from the best fit of the SSPM solutions to 

the exact Mie results in spherical case at wavelengths 8.5, 11.0, and 12.0 pm. 
The values of SSPM extinction efficiency corresponding to the best fitting A a e 
are also shown in comparison with the Mie results. The parameter A a e are 
presented in an absolute scale (in the unit of incident wavelength) and also the 
percentage with respect to the particle radii. 

Fig.8 As the same as Fig.7 except for absorption efficiency. Pronounced tunneling 
effect is evident from the result for 12.0 pm wavelength, as the absorption 
efficiency is substantially larger than unity at particle radii near 10 pm. 

Fig.9 The absorption of hexagonal ice crystals computed from various methods. The 
particle sizes are defined by their maximum dimensions with an implicit 
inclusion of aspect ratio defined in Eq. (1). The expansion parameter A a e and 

A a a required for SSPM calculation in the hexagonal case are approximated by 
their counterparts in a spherical case, that is, the number of layers of stretched 
potential and their thickness in hexagonal case is the same as in the spherical 
case. 

Fig. 10 As the same as Fig.9 except for extinction efficiency. 

Fig.l 1 The combination of the FDTD and refined SSPM results for the absorption and 
extinction efficiency for the particle sizes ranged from 1 to 10000 pm. A 
smooth transition between the two methods are evident. 
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Fig. 12 The effective sizes and mean maximum dimensions computed for the 30 size 
distributions that are taken from Fu [1996] and Mitchell et al. [19996]. 

Fig. 13 The variation of the volume-normalized extinction coefficient versus the 
effective size for 12 wavelengths in the IR window region. 

Fig. 14 The mean extinction efficiency calculated for the 30 size distributions (circle 
symbols). The solid lines are the parameteization results obtained from fitting 
the “exact” data in terms of the second order polynomial function of 1 / D e . 

Fig. 15 As the same as Fig. 14 except for mean single-scattering albedo. An increasing 
trend is evident for the variation of the single-scattering albedo versus 1 / D e at 
wavelengths 8.0, 8.5, 9.0, 9.5 and 10.0 pm. The opposite is noted for the other 
wavelengths. 

Fig. 16 The variation of the asymmetry factor versus D E for the 12 wavelengths (circle 
symbols). The solid lines are the parameterization results based on power law 
fitting 
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Table 1. The fitting coefficients for the parameterization of the mean extinction 


efficiency, which is defined by the Eq.(22). 


Wavelength (pm) 

T\o 

Til 

^2 

8.0 

1.989 

4.336 

-2.776E+1* 

8.5 

1.979 

6.117 

-4.877E+1 

9.0 

1.982 

5.219 

-4.248E+1 

9.5 

1.992 

3.891 

-3.743E+1 

10.0 

1.984 

5.208 

-8.194E+1 

10.5 

2.032 

-3.643 

-4.075E+1 

10.8 

2.036 

- 6.777 

5.749 

11.0 

2.029 

-5.179 

3.789 

11.5 

2.012 

-1.431 

-1.061E+1 

12.0 

2.006 

1.283 

-2.800E+1 

12.5 

2.000 

2.787 

-3.323E+1 

13.0 

1.998 

3.201 

-2.664E+1 


*: 2.776E+1 indicates 2.776 xlO +1 . 
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Table 2. The fitting coefficients for the parameterization of the mean single scattering 


albedo, which is defined by the Eq.(24). 


Wavelength (pm) 

C„ 

Cl 

C 2 

8.0 

5.294E-1 

2.536E-1 

1.541E+1 

8.5 

5.205E-1 

1.421 

7.036 

9.0 

5.198E-1 

1.178 

9.676 

9.5 

5.191 E- 1 

1.040 

1.063E+1 

10.0 

5.116E-1 

1.770 

-1.811 

10.5 

5.230E-1 

-8.619E-1 

1.271 

10.8 

5.384E-1 

-2.454 

1.225E+1 

11.0 

5.475E-1 

-2.552 

1.417E+1 

11.5 

5.562E-1 

-1.974 

9.011 

12.0 

5.633E-1 

-1.919 

9.232 

12.5 

5.556E-1 

-1.605 

7.246 

13.0 

5.678E-1 

-2.051 

1.237E+1 
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Table 3. The fitting coefficients for the parameterization of the mean asymmetry 


parameter, which is defined by the Eq.(25). 


Wavelength (pm) 

<; 

K 

8.0 

8.317E-1 

3.279E-2 

8.5 

8.217E-1 

3.565E-2 

9.0 

8.399E-1 

3.120E-2 

9.5 

8.560E-1 

2.779E-2 

10.0 

8.917E-1 

2.092E-2 

10.5 

9.27 IE-1 

1.387E-2 

10.8 

9.092E-1 

1.771E-2 

11.0 

9.080E-1 

1.501E-2 

11.5 

8.922E-1 

1.394E-2 

12.0 

8.519E-1 

2.355E-2 

12.5 

8.553E-1 

1.932E-2 

13.0 

8.290E-1 

2.682E-2 
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